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Abstract 

We present together in this document the paper B. Garcia- Archilla, 
Shiskin mesh simulation: A new stabilization technique for convection- 
diffusion problems, Comput. Methods Appl. Mech. Engrg., 256 (2013), 1- 
16, and the manuscript B. Garcia- Archilla, Shiskin mesh simulation: Com- 
plementary experiments, which contains some of the experiments that, for 
the sake of brevity were not included in the first one. Thts second paper 
with complementary experiments is neither self-contained nor intended to 
be published m any major journal, but is intended to be accessible to any- 
one wishing to learn more on the performance of the SMS method. Follow- 
ing a principle of reproducible computational research, the source codes of 
the experiments in the present papers are available from the author on re- 
quest or on "http: //personal .us . es/bga/bga_files/sof tware_bga.html"' . 
Excluded are the codes corresponding to Example 7 which is subject of 
current resarch 
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Shishkin mesh simulation: A new stabilization 
technique for convection-diffusion problems 

Bosco Garcia- ArchillaQ 



Abstract 

A new stabilization procedure is presented. It is based on a simulation 
of the interaction between the coarse and fine parts of a Shishkin mesh, 
but can be applied on coarse and irregular meshes and on domains with 
nontrivial geometries. The technique, which does not require adjusting 
any parameter, can be applied to different stabilized and non stabilized 
methods. Numerical experiments show it to obtain oscillation-free ap- 
proximations on problems with boundary and internal layers, on uniform 
and nonuniform meshes and on domains with curved boundaries. 

Key words. Convection-dominated problems. Stabilized methods. 
Finite-element methods. Galerkin method. SUPG method. 

1 Introduction 

The numerical solution of convection-diffusion problems when convection dom- 
inates is, despite more than 30 years of research, a challenging problem nowa- 
days. Standard finite-element or finite-difference methods typically suffer from 
unphysical or spurious oscillations unless meshes are taken so fine that are use- 
less for all practical purposes. The reason is the presence of layers or thin regions 
where solutions change fast. Modification of standard methods, known as sta- 
bilized methods have been proposed in the literature, from upwind methods 
35 years ago [45], to strongly-consistent stabilized methods like the stream- 
line upwind/Petrov-Galerkin (SUPG) method [7], also known as the stream- 
line diffusion finite element method (SDFEM), or the Galerkin least-squares 
(GALS) method [21]. More recently, local projection stabilization (LPS) meth- 
ods, [1], [B], HH], continuous interior penalty (CIP) methods 0, or discontinuous 
Galerkin (DG) methods [23], [35] have been introduced, to cite a few of the many 
techniques proposed (see [36], [38] for a survey of methods). It must be noticed, 
however, that computational studies (see e.g., |2], [26]) find it hard to put a par- 
ticular method above the others. It must be also mentioned that most of these 
methods depend on at least one parameter about which there is no unanimous 
agreement on its optimal choice in practical problems |27| . 
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A different approach is to use layer-adapted meshes. Among these we cite 
Shishkin meshes (described below) [31], [37], which have received considerable 
attention in recent years [TJ, [H], [H], [15], [SI], [33] [32], [13] and [1?]. How- 
ever, it is generally acknowledged that the main drawback of Shishkin meshes 
is the difficulty to design them on domains with nontrivial geometries, although 
some works overcoming this difficulty can be found in the literature [35] , [55] • 

The method we propose, however, does not suffer from the above indicated 
drawbacks: It does not depend on parameters and, although it is based on the 
idea of simulating a Shishkin mesh, the experiments we present show it produces 
excellent results on domains with nontrivial geometries. 

We consider the problem 

- eAu + b-\/u + cu = f, in n, (1) 
u = gi, vad^D, ^ = .g2, in 917 at. (2) 

Here, $7 is a bounded domain in R'^, d = 1,2,3, its boundary dO, being the 
disjoint union of T £> and F^r, 6 and c are given functions and e > is a constant 
diffusion coefficient. We assume that T~ C dVlu^ T~ being the inflow boundary 
of C i.e., the set of points x € dU, such that h{x) ■ n(x) < 0. 

It is well-known if e <C sup{|6(a;)| | x e f2} (|-| being the euclidean norm) 
boundary layers are likely to develop along (9f7\r^, although they have different 
structure on r" = {x e dD, \ b{x)-n{x) = 0} and r+ ^ {x e dQ \ b{x)-n{x) > 0}. 
As already mentioned, these boundary layers, when present, are responsible of 
the spurious oscillations that pollute the numerical approximations obtained 
with standard methods unless extremely fine meshes are used. For uniform 
meshes, oscillations typically disappear when the mesh Peclet number 

||&|Il°°(o)2 h 
2e 

{h being the mesh size) is of the order of 1. 

Let us briefly describe now the idea of the method we propose in the following 
simple problem: 

L{u) = -eu"{x)+b{x)u'{x)+c{x)=f{x), < a; < 1, (3) 
m(0) = m(1) = 0. (4) 

In ([3]) we assume that b, c and / are sufficiently smooth functions, and that 

< B < min b(x), < min c(x). (5) 

xe[o,i] xelQ.i] 

The standard Galerkin linear finite-element method for ([S]!!]) on a partition or 
mesh = xo < xi < . . . < x,7 = 1 of [0, 1] obtains a continuous piccewise linear 
approximation Uh{x) to u. As it is customary, h denotes the mesh diameter, 
h = maxKjxj hj, where hj — Xj — a^j-i, for j ~ 1, . . . , J. The approximation 
can be expressed as U{x) = uiipi{x) + • • • + uj-iipj-i{x), where the (pj{x) are 
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the basis or hat (piecewise hnear) functions taking value 1 at the node Xj and 
in the rest of the nodes of the partition (thus, U(xj) = Uj). The values Uj, 
j — 1, . . . , J — 1^ are obtained by solving the linear system of equations 

a{uh,ip,) ^ {f,ipi), i = l,...,J-l, (6) 

where, a is the bilinear form associated with which is given by 

a{v, w) = e(v' , w') + {bv' + c, w), 

(•, •) being the standard inner product in L^(0, 1), 

[La) = I f{x)g{x)dx. 
Jo 

The Shishkin mesh with J = 2N nodes is composed of two uniform meshes 
with N subintervals on each side of the transition point xn — 1 — a, where 

1 2 

<T = min(-, -e log TV), 

for an adequate constant /3, that is, Xj = j{l — (t)/N, for j — 0,. . . ,N, and 
XN+j — xn -\- jcr/-^, for j — 1, . . . , iV. Let us consider the coarse and fine grid 
parts of the Galerkin approximation given by 

Uc{x) ^ UiLpi{x) H h UN-l(pN~l{x), 

Uf(x) = Un+WN+i{x) H 'rU2N-W2N^l{x), 

SO that Uc + UNifN + Uf is the Galerkin approximation on the Shishkin mesh. 
Since for i = 1, . . . , J — 1, the support of the basis function ipi is [xi^i, Xi^i], 
we have a(Uc, fN+j) — and a(Uf, (pj) = 0, for j = 1, . . . , — 1. Consequently 



the system ^ on the Shishkin mesh can be rewritten as 

a{Uc,V^) ^{f,Vi). ^ = l,...,A^-2, (7) 

a{Uc.,fi)+UNa{LpN-,^i) ^{f,(pi), (8) 

a{Uc, ipn) + UNa{(pN,ipN) + a{Uf, (Pn) = (/, Pn), (9) 



UNa{(pN,Pi) + a{Uf,pi) ^{f,(pi), i = N + 1, . . . ,2N - 1. 

(10) 

We notice that were it not for the presence of the UNai^N tPi) in (El), the 
system ([THH]) would be the equations 

a(U,p,) = {f,p,), z = l,...,iV-l, (11) 

of the Galerkin approximation U = Uiipi + • • • + Un-iP>n-i for the problem 

~eu"{x)+b{x)u'{x)+c{x)= f{x), 0<x<l-a, (12) 
u{0) = u{l ~a)=0. (13) 
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Figure 1: Galerkin approximation on a uniform mesh with iV = 9 (continuous 
line) to the solution of with e = 10'^, a = 4elog(2 * N), b = f = 1, 

and c — 0. The Uc part of the Galerkin approximation on a Shishkin mesh with 
J = 18 (broken line) for same e and /. Circles are the values of the true solution 
on the nodes. 

The Galerkin approximation U for this problem, unless sN > 1/2, is likely to 
have spurious oscillations of large amplitude as we show in Fig. [l]for e = 10~^, 
a = 4elog(J), b{x) — f{x) — 1, c — and = 9. It is however the presence 
of uno.{'Pn , ^i) in. equation that suppresses the oscillations, as we can see 
in Fig. [1] where the component Uc of the Galerkin approximation on a Shishkin 
grid with J = 2N = 18 is also shown (discontinuous line) together with the true 
solution at the nodes of the coarse part of the mesh. 
It is remarkable that just by adding the value 



to the last equation of the Galerkin method for ([T2HT3]) we get the oscillation- 
free approximation Uc- Obviously, in order to have the value of a* we have 
to solve the whole system (fZl fTUl) . In the present paper, we introduce a tech- 
nique to approximate a* without the need to compute the whole approximation 
on the Shishkin grid. In Fig. [TJ the approximation computed with the esti- 
mated a* is indistinguishable from Uc- Numerical experiments in the present 
paper show that, in two-dimensional problems, the oscillation-free approxima- 
tion on a coarse mesh can be obtained by this technique at half the computa- 
tional cost of a Shishkin grid, and a more substantial gain can be expected in 
three-dimensional problems. 

Furthermore, this technique can be extended when the grid is no part of 
any Shishkin grid, while, at the same time, managing to get rid of the spurious 
oscillations. This allows to obtain accurate approximations on domains with 
non trivial boundaries, where Shishkin meshes may be difficult to construct. In 
spite of this, we call the new technique Shishkin mesh simulation (SMS), since 
it was derived, as described above, in an attempt to simulate Shishkin grids. 

We must mention, however, that in the present paper we only consider the 
case of dominant convection, both in the analysis and in the numerical experi- 



a* = UNa{^pN-i,'PN) 



(14) 
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ments. The question of how to modify the method (if necessary) when the mesh 
Pcclet number Pe tends to one wih be addressed elsewhere. 

It is wcU-known that the Galerkin method is a far from ideal method in 
convection-diffusion problems. Let us also notice that despite the good proper- 
ties of stabilized methods developed in recent years, the SUPG method is still 
considered the standard approach . For this reason, in the numerical exper- 
iments presented below, we compare the new method with the SUPG method. 

The rest of the paper is as follows. In Section[2]we describe the SMS method 
in detail. In Section [3] we present a limited analysis of the new technique. 
Section 21 contains the numerical experiments and Section [S] the conclusions. 

2 The new technique: Shishkin mesh simulation 
2.1 The one-dimensional case 

We consider ([3H11) satisfying ([S]). Given a partition ^ xq < xi < . . . < xj ~ 1, 
of [0, 1], we denote by Xh the space of continuous piecewise linear functions, and 
by Vh the subspace of Xh of functions taking zero values at x = and x = 1, 

so that we can express 

Xh = span{(y5o} ® V/i © span{v3j}. 

We consider the operator Lh given by 

Lh{vh) = bv'h + cvh, Vh e Vh- 

(See also Remark [T] below). We denote by Uh G Vh the standard Galerkin linear 
finite-element approximation, and by Uh the SMS approximation. This is found 
as the solution of the least-squares problem 

- 3^"^ ^J^h{uh) ~ fWmo.xj^,) ^ (15) 

subject to the restriction 

a(u/,.,<P/0 + "'Pft(2;j-i) = (/, <^), <p/i e T4. (16) 

Observe that since kphi^^j-x) — for kph not proportional to i^aj-i, and recalling 
how we defined a* in (|14l) , equation (|T6)) is similar to ([THH]) , and then, the least- 
squares problem ()15p is the way of finding the value a hopefully close to a* . 

Notice also that the restriction (fT6|) is, in fact, a set of as many independent 
restrictions as interior nodes or nodal basis functions in Vh- Consequently, in 
the optimality conditions, there must be a Lagrange multiplier for every interior 
node. We gather all this multipliers in a function Zh € Vh whose value at every 
node is that of the corresponding Lagrange multiplier. Thus, the optimality 
conditions of the SMS approximation can be written as the following linear 
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problem: find iih, E Vh and a G R such that 



Zh{x,J-i) 



(17) 

= 0, (18) 

= (/,^), cpheVh, (19) 



where here and in the sequel denotes the standard inner product 

in L^I). 

Remark 1 Observe that for linear elements the operator coincides with 

J 

Lhivh) = ^(-e< + bv'^ + CT/,)|,^^_^,^^, , (20) 

that is, L applied element by element. This expression of is better suited 
to the SMS method for higher-order elements, a topic that will be studied else- 
where. 

Remark 2 The fact that the approximation Uh is found by solving the least- 
squares problem ([15]) may suggest a possible relation with the GALS method 
(compare for example least-squares problem ([T5|) with that in [3BJ p. 327]). 
However, they are very different methods, since for the examples in Section [?] 
in the present paper, SUPG and GALS methods are identical (see e.g. [21]), 
and, as shown in Section H] the SMS method and the SUPG method produce 
markedly different results. 

2.2 The multidimensional case 

We will assume that b{x) ^ for x G il, and that every characteristic (i.e., 
solution of dx/dt = b{x)) in enters and leaves in finite time. Also for sim- 
plicity, we will assume that the domain in ([T]^ has a polygonal or polyhedral 
boundary. 

Let Th a triangulation of it, that is, a partition of the closure of in 
n-simplices with only faces and vertices in common. For every t S 7h, let A/'(t) 
be the set of its vertices, and let J\fh — Ur&Th-^i''') be the set of vertices of Th- 

Similarly to the one-dimensional case, let Xh the space of continuous piece- 
wise linear polynomials. We express 

Xh=X-®Vh(SX+, 

where (ph = on Fd if iph G Vh, and for iph in (resp. Xj^), if (ph{x) 7^ 
for X £ Afh, then x G F^ (resp. F£)\F^). In the standard Galerkin method, 
first an element G X^^ ® X'^ is selected such that the restriction uj^ to 
Fu is a good approximation to the Dirichlet data 51 in ([5]). This restriction is 
typically the interpolant or the L^(F£))-projection onto the restrictions to F^i 
of functions in X/^ . Then, the Galerkin approximation u/j G uf^ + satisfies 

a{uh,(ph)^{f,'Ph)+e{g2,^h)r,^, V(p/i G Vh, (21) 
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where here and in the sequel, (•, •)r denotes the inner product on F C dft, 
and 

a{v, w) — e{\7v, \/w) + [b ■ Vu + cw, w), v,w £ H^{D,). 
For the new method, similarly to the one-dimensional case, we consider 

Lhivh) = b ■ Vvh + cvh- 

In order to describe the new approximation, we must set up the multidimen- 
sional version of the last interval (x,/_i,a;j) in the one-dimensional case. For 
this purpose we denote 

= (r+uF")nFr,. 

We consider a suitable with F^ C d^l^ (to be specified below) and let us de- 
note by Ns and the set of vertices in drij^\T]j and their indices, respectively, 
that is 

Ms = N{dn+\TD), = {j e N I x, e Ms}, 

and by R'*'"' we will refer to the set of real vectors of the form (tj)jgN^ . The points 
in Ms and the set $7^ will play a role similar to that of xj-i and {xj-i^xj) in 
the one-dimensional case. This can be seen in Fig. [2l where, shadowed in grey, 
we show the set fl^ for a triangulation of the unit square for b = [1,1]"^, with 
Tu = Sfi, so that r° U r+ is consists of the sides y = I and a; = 1. We show 
the points in A/^ marked with circles. 




Figure 2: A triangulation of the unit square with F'^UF+ marked with a thicker 
line, the set in grey and the the points of Ns with circles. 

The approximation Uh G + Vh, is then found by solving the least-squares 
problem 

min ||L,i({t/i) - /II 2,„, (22) 

subject to the restriction 

a{uh,iph) + ^ tj(ph{xj) ^ {f,iph) + e{g2,iph)r„, fh ^ Vh- (23) 
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The optimality conditions of this problem can be written as the following linear 
problem: find Uh E + Vh, Zh G Vh and t G R^* such that 

{Lh{uh),Lh{iph)) ^2^^'y^+-,- ahifh,zh) = {Lh{^h)J)L2(n\n+y (Ph e Vh, (24) 

zh{x,)^0, jeNs, (25) 

a{uh,(ph) + ^ tj(ph{xj) ^ {f,(ph) +e(52,<P/i)rjv) 'Ph e Vh, 

(26) 

where, as in the one-dimensional case, z^ is the Lagrange multiplier of restric- 
tion ([23]) . 

Let us specify now the set J . The obvious choice of setting = where 

Bk = U ^' 

rn(r^+)7^0 

as in Fig. [21 may lead to an unstable method if there are nodes x S Afh interior 
to Bfi, as it can be easily checked in a numerical experiment. To understand why 




Figure 3: A triangulation of the unit square with F'^ U F+ marked with a thicker 
line and the points of Afs with circles. In grey, the sets Bh (left plot) and fi^ 
(right plot) for b = [1, 1]^. 

in such a case the SMS method may be unstable, consider the limit case e = 0, b 
constant and c = (so that the bilinear form a is skew-symmetric) and consider 
a mesh as that depicted in Fig. [3J where there is one node Xi interior to Bh 
who has only one neighbour node Xk G Afs- Then, on the one hand, the basis 
function (fi of the node Xi satisfies that ah{ipi,(ph) — for all (ph E Vh, except 
when ifh is a multiple of the basis function (pk of the node Xk- Furthermore, 
since the support of ipi is contained in Bh, we have that \\Lh{<p>i)\\L2(^Q\g^-j = 0. 
Consequently, taking Vh — Pi, Zh — 0, ti — for I ^ k and tk = —a{ipi, ipk) we 
have a nontrivial solution of 

[Lh{vh),Lh{'Ph)) ^2^^^\f^-i+)- a,h{ph,Zh) iph&Vh, (27) 

Zh{xj)=0, jeNs, (28) 
a{vh,(ph) + ^ tjiph{xj) =0, iphE Vh, (29) 
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if we set = Bh. It can be easily checked that this is also the case if the 
node Xi interior to Bh is connected to more that one node in A/^. For e > 0, one 
can expect a unique solution but, as we have checked in practice, it is a solution 
where u{xi) = 0{l/e). 

To avoid this situation we remove from B^ the upwind triangle of any interior 
point. More precisely, for every node Xi G Afh, let us denote by r_(xi) its upwind 
triangle, that is 

T^{x,) = {TeTh\{x^-Xb}nT^^!), A^O}. 

Then we define 

o 

If there are two upwind triangles because Xj — Xb is and edge for A smaller than 
some Aq, we may select the first one of them after ordering all the elements. 

With this choice of il'^ , as we will show in Section 13.2.11 the few cases were 
there are nontrivial solutions of ((27l - [29| have an easy fix from the computational 
point of view. We also note that a slightly different choice of is taken 
in Example [7] in Section HI where the vector field b has closed integral curves. 

2.3 Further Extensions 

In the same way that the Shishkin meshes are not restricted to the standard 
Galerkin method, the new technique, originally motivated by Shishkin meshes, 
is not restricted to the Galerkin method either. In this section we comment on 
how to use it with some other methods of widespread use, such as the SUPG, 
GALS, LPS and GIF methods. 

In these methods the approximation li/, G u^ + V/j, instead of satisfying (j2ip . 
satisfies 

ah{uh,iph) ^ {f,Vh)h + £{92,^h)ri^, y^h^Vh, (30) 

where and (•, ■)h are mesh-dependent bilinear forms. For the SUPG method 
they are given by 

ah{vh,^h) = a{vh,iph) + ^ Sr{L\^{vh),b ■ iph)T, (31) 

reTh 

if, '^h)h = (/, <y5ft) + X! ^ ' (^^) 

where 5r is and adequately chosen parameter and denotes L restricted to r. 
In the GALS method the term b ■ Vfu is replaced by L\^{tfh) (see e. g., [33] 
for a full description of these methods). The extension of the SMS technique 
to these methods is as follows: we solve the least-squares problem (|22|) but we 
replace the restriction ((23)) by 

ah{uh,<^h)+ ^ tjtphixj) = {f,iph)h + £{g2,Vh)rN- G Vh- (33) 
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In the numerical experiments in the present paper, we only consider the SMS 
for the SUPG method. To distinguish it from the SMS method in the previous 
section, we will refer to them as Galerkin-based and SUPG-based SMS methods. 

We could also consider the extension to higher-order finite-element methods. 
The extension seems straightforward since it consists in replacing the approxi- 
mation space Vh by piecewise quadratics or pieceswise cubics (together with 
in ((20)) '). Numerical experiments (not reported here) indicate that this straight- 
forward extension does not give as good results as the piecewise linear elements, 
at least for one-dimensional problems. The reason seems to be the need to 
redefine the set fi^ for higher-order elements. This will be subject of future 
studies. 

3 Analysis of the SMS method 

3.1 The one-dimensional case 

Some understanding of why the new method is able to suppress or, at least, 
dramatically reduce the spurious oscillations of standard approximations can 
be gained by analyzing the problem ([SH!]) when 6 is a positive constant and 
c = 0. We do this first for the Galerkin method. 

The Galerkin method. 

We first consider the limit case e = 0. In this case, it is easy to check that 
the Galerkin approximation satisfies the equations 

b 

-{uj+i - Uj-i) ^ fj = f{x)ipj{x)dx, j = l,...,J-l. (34) 

Thus, summing separately odd and even-numbered equations we obtain the 
following expressions, 

2 ^ 

=uo + -^/2.-i, J = !,...,( J' -l)/2, (35) 

2 — 1 

2 (J'-i)/2 

U2j-i=uj'^- J2 j = l,...,(J'-l)/2, (36) 

where, here and in the rest of this section 

^, f J, if J is odd, 

1 J — 1, if J is even. 

We notice that when J is odd, the expressions in (jHSHMl) are the discrete 
counterparts of the problems 

bu^ = f, < a; < 1, u(0) = 0, (37) 
bu^^ f, < a; < 1, u{l) = 0, (38) 
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which, unless / has zero mean, have different solutions. Thus, since for suffi- 
ciently smooth / the expressions ((35H36|) are consistent with ([37H38| , oscillations 
in the numerical approximation are bound to occur whenever / is not of zero 
mean. For J even, it is easy to check that the Galerkin equations have no 
solution unless /i + /a + • • • + fj-i = 0, in which case the solution is not unique. 
Thus for J even the Galerkin method is not stable. 

However, for J odd, the Galerkin method is stable in the following sense 



where 
and 



\u,\\^<^\\f\\_^, (39) 



-h^WShW^, (40) 



j (./'-l)/2 
52,=^/2.-l, ^2,-1= J = l,...,(J'-l)/2. (41) 

The SMS technique. Here the equations are 

^{uj+i-Uj^i) + aSj^i{j) = fj j = l,...,J-l, (42) 
where 6j-i denotes Dirac's delta function 

^■'-^(■^)-[ 1, ifj = J-i. 

As Proposition [T] below shows, the SMS method is stable independently of 
the parity of the grid. To prove stability we will consider the auxiliary function 
qii G Vh whose nodal values are 

1 - f-l)-'' 

q,= j = l,...,J-l. (43) 

It is easy to check that when J is odd qh satisfies 

a{qh,iph) = ^h{x.j-i), ifh e Vh, (44) 
whereas when J is even it satisfies 

a{qh,fh) = 0, (fh. e Vh- (45) 

Proposition 1 There exist a positive constant C such that for e = the SMS 

approximation satisfies 

\\uH\\^<U\\f\\-H + ^Ar\\, (46) 



6J 



where r = {f, Lhiqh))mo.xj^i)- 
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Proof. Recall that the SMS approximation Uh is part of the solution (u^, a, Zh) 
of system p7HT9)) . Due to (gS]) and ^ we have that a{qh,Zh) 0, so 

that taking iph ~ qh in (|17p we have 

{Lh{uh), Lh{qh)) ^^f^ ,, ^^^j if,Lhiqh))L^o,xj^^)- (47) 

On the other hand, we notice that L{qii) — 2{—iy/hj, for x e {xj^i,Xj) and 
j = 1, . . . , J — 1, so that equation (|47|) becomes 

2g(_i).5^^-^ = 2X:(-l)^7^ r /(x)rf^- (48) 

j=l "■■7 j=l '^J 

To prove we treat the cases J odd and J even separately. For J odd, we 
notice that since the equations in are those of the Galerkin method except 
for the last one, applying the stability estimate ([551) we have 

||^i/.IL<^(||/IL, + |a|). (49) 
Also, since Uh = Uh — aqii, from (|48|) it follows that 



r being the right-hand side of (|48j). Thus, applying again the stability esti- 
mate (p9|) . we have 

|a|<2||/|L, + ^^^M. (50) 

From this inequality, and the fact that 2(J — 1) > J, the stability esti- 
mate (|46)) follows. 

For J even, taking (^/i = in (fT9|) and applying (j45|) we have 



2 

= {b{uh)x - /, g/i) + agj-i = aqj-i + ^ X! /2j-i, 

so that 

J/2 

a = E/2.-i- (51) 



Also, since the equations in p2|) are those of the Galerkin method except for 
the last one, which does not appear in ([55]) and ([551) for J even, we have that 
for J even ([55)) and also hold with m/j replaced by u/j. Thus, it follows that 

\\uh\\<l\\.f\U + \uj-i\. (52) 
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The value of uj-i is then obtained from (|48p by replacing Uj by their values 
as expressed in ([35H36|) . which gives 



and, thus, 

|u.-i|<^ll/IL. + ^^H, (53) 
which, together with (|52]) shows that (|46| also holds for J even. □ 



Once the stability is proved, we investigate the convergence of the SMS 
method by considering the basic solution, 

bu^ = /, u{0) = 0, 

and the error 

eh = Uh - Ih{u), 

where Ih{u) is the interpolant of u in Xh. We prove different estimates depend- 
ing on whether 

max - = Cft.^, (54) 

holds for some constant C > 0. 

Proposition 2 There exist a positive constant C such that for e = the error 
= Uh — Ih (u) of SMS approximation satisfies the following estimates: 

||g.L<C/i||/'||^,(„_i), (55) 

and, if |5-^[ j holds, 

||e~.L<C;i2(|l/'IL + ||/"ILi(o,i))- (56) 

Proof. We have 

b,^ ^ . / b 



-(gj+i -ej_i) + (^a- -u(l))(5j_i(j) =rj, j = 1,...,J -1, (57) 



where 



Xj-l 

Since ^ 

/ ' iblhiu)^ -f)dx = 0, j = 1, . . . , J - 1, 

so that {b{qh)xMh{u)x - f)L^(o.xj-i) = 0, it follows that 



ib{^h)x,b{eh)x)mo.xj^i) =0. (58) 
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Thus, applying (^5)) with / replaced by r = [ri, . . . , tj_i] and r by 0, we have 

\\eH\\^<l\\rH\\_f,. (59) 



Integration by parts reveals 



n^){f{\-Vj)dv) 



dx. 



and, lurther, 



dx. 



Then, it is easy to check that 

llr||_,<C/.||/'||^,(o^i), (60) 



or 



Irll , < max 

fe<(J-l)/2 



J=fe 
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+ C•/^^(ll/'ILoo(o,l) + ll/"llLMoa))• (61) 

We notice that from ([5^ and ([50]) the estimate ([55]) follows. On the other hand, 
if (IMl) holds, then the right-hand side of §^ is 0{h'^), so that ^ and (HH) 
imply (HH). □ 

Let us finally comment the case e > 0. In this case, equations (j42|) should 
be modified by adding the term 

hj+i hj 

to the left hand side. Then, the stability estimate holds but with 
on the right-hand side replaced by 

i<j<.j Hj 

added to the right-hand side. Thus, we have the same stability estimate (|46p 
with the factor 6/6 replaced by 12/6 whenever 



6 

48 J i'<j<j' 



Finally, with respect to the error, standard asymptotic analysis (see e. g. [551 
Theorem 1.4]) shows that the solution u of ([3][4]) and the basic solution uq 
(solution of (1571) ') satisfy that 

max \u(x) - uo(x) + wo(l)e^^'=~^^/1 < Ce, 



14 



for some C > independent of e. Since e''*^^ ^'^^^ < e for a: < 1 
then, whenever 

, ^ £|log(g)l 

the bounds (l55|) and hold if ((Ml) also holds. 

Remark 3 As mentioned in the Introduction, in the present paper we are con- 
cerned with large Peclet numbers. For this reason and for simplicity, we have 
not pursued conditions which allow for larger values of e than (15^ . Obviously, 
a more detailed (and lengthy) analysis would allow to obtain more general con- 
ditions on e. 



-£|log(e)|/6, 
(63) 



3.2 The multidimensional case 
3.2.1 Uniqueness of solutions 

In this section we limit ourselves to study the possibility loosing uniqueness in 
the new method when e = and how to overcome this in practice. These cases 
are an indication of when to expect the method to perform poorly for small 
e > 0. In the rest of the paper, we denote 

We start by characterizing the solutions of (|27l - [29l) when e = 0. As we now 
show, they are given by those elements vu G Vu satisfying 

II^'^(«'OIIl.(o,)-0. (64) 
To see this, we take iph = vu in ([27| and get 

||i/i(w/i)|i^2(o^) +a{ih,Zh) = 0. 
But taking iph — Zh in (1291) we have 

a{vh,Zh) = - ^ tjZhixj) = 0, 

where in the last equality we have applied p8|) . Thus, ([64)) follows. Conversely, 
if (|64l) holds, this implies Lh{vh) = on But since we are assuming that 
£ = 0, we also have that a{vh,Wh) = {Lh{vh),'Wh) for any VhjWh S Vh- Thus, 
if for a basis function ipj we have ^ a{vh,yyj) — {Lfi{vh), fj), it follows that 
its support cannot be entirely inside fi/j, so that Xj must be in Ms- Then 
taking Zh = and tj = —a{vh, (pj) we have a solution of (|27H29| . 

We devote the rest of this section to comment on a case that may easily arise 
in practice where there is a Vh not entirely null in ft that satisfies ([M)) (and, 
hence, the lost of uniqueness in the SMS method when e = 0). We will also 
comment on what possible remedies can be applied in practice in these cases. 
In the rest of the section we assume c — 0, that is, Lh{wh) — b ■ 
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Figure 4: Two triangulations of the unit square with of disconnected interior. 
The set fi^ for 6 = [1, 1]-'" is shadowed in grey and points oijVs are marked with 
circles. 

Consider an example like those in Fig. |4l that is, when the interior set flh 
of Qh is not connected. We may have Vh G Vh not null in fl and yet satisfy- 
ing (|64p . To see why, take for example the left plot in Fig. |4]and assume that 
b is constant with strictly positive components for simplicity. Let Xj and Xk 
be the most downwind vertices of the isolated triangle t of Clh. Observe that 
since their associated basis functions tpj and ipk have linearly independent gra- 
dients, it is possible to find real values a and /3 so that a^pj -\- f^ipk 7^ but 
b ■ {aS/ipj + l3Vipk) — 0. Setting then Vh = aipj -f Pip^ we have that Vh ^ Q 
but it satisfies (|64|) . If b is not constant, then it is possible to find Vh such that 
II^'iIIl2(o) = 1 but \\Lh{ 

"^i^^W L^{hh) ^ diam(r) ; there may not be nontrivial 
solutions of (|27l - [29l) but the method may not be very accurate since residuals 
of size diam(r) allow for errors of size 1 = ||^/i|Il2(q') . Similar arguments show 
that this is also the case of the right plot in Fig. |4l These arguments are easily 
extended to tetrahedra in three-dimensional domains. In Section [3.2.21 we show 
some other cases where there are nontrivial Vh satisfying We also state 

general conditions to prevent them. 

From the practical point of view, all those conditions are satisfied if the 
meshes are designed with a strip of elements on as indicated in Fig. [5l One 




ponents in Vlh- The set 57^ for b = [1,1]"^ is shadowed in grey. 



possible way to form this strip is to take the nodes on the outflow boundary 
(which we assume partitioned into segments or faces) and displace them along 
the normal vector to obtain the interior nodes of the strip (see Fig. Then, 
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sides or faces in the outflow boundary are replicated in the displaced nodes, so 
that rectangles are obtained in two-dimensional problems and triangular prisms 
in three-dimensional ones. In the two-dimensional case, the triangulation of the 
strip is obtained by joining two opposite vertices of the rectangles. In three 



Figure 6: The process to build a strip of elements along T j^ ■ 

dimensions, tetrahedra can be obtained for example by adapting to triangular 
prisms the so-called Kuhn's triangulation of the cube, [3], [T7], [55]. This is a 
partition of the unit cube into 6 tetrahedra, all of them having the origin and 
the vertex of coordinates (1, 1, 1), as common vertices. It is a simple exercise to 
adapt Kuhn's triangulation of the cube to triangular prisms. 

We notice, however, that there are times in practice when one cannot choose 
the grid because it is part of a larger problem, or because it is designed to 
satisfy some other constraints, etc. In Section fS. 2. 2) we comment on a different 
procedure to avoid isolated components in flh- 

3.2.2 Further cases of lack of uniqueness 

We complete the study of the previous section on the cases where there are 
nontrivial v/^ G Vh satisfying ([M)) . We assume £ = and c — Q here as well. We 

start by checking that only Vh — Q satisfies ([M)) when Clj^ is disconnected due to 
an element upwind of a node interior to Bh, as it was the case discussed in Fig.|21 
Let us take that example again. We first notice that since we are assuming that 
e = 0, then, ((64)) imply that vt is constant along characteristics, and since 
Vh = Q on r_, it must vanish on the big component of Clh- This implies that 
Vh also vanish on the most upwind vertex xu of the isolated triangle r. It also 
vanish on the vertex on the boundary. So it is only on the remaining vertex Xi 
where it may not vanish. But recall that the triangle r is upwind of Xi , so that 
h ■ W(pi 7^ 0. Thus, Lfi{vh) = on r implies that Vh — on t. 

However, even when Clh has a connected interior, there may be nontrivial 
Vh G Vh satisfying when there are triangles (resp. tetrahedra) in Clh, down- 
wind of and with one side aligned with the (constant) wind velocity b (resp. 
one face parallel to b). This is the case of the triangles with one side plotted 
with a thicker line in the grids in Fig. [7] for b = [1,1]"^. The basis function of the 
vertices opposite those sides have their gradients orthogonal to b. For example, 
for the triangle t marked with a + sign in the center plot of Fig. [71 let Xi be the 
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Figure 7: Triangulations of the unit square with triangles in Clh downwind of il^ 
and with one side (thicker hne) parahel to the wind velocity b = [1,1]"^. The 
set fi^ is shadowed in grey. 

vertex opposite to the shadowed side on the triangle, and Xj and Xk the other 
two vertices. Since b ■ V (pi = in r, and (^j = on the rest of fi/j, then Vh = ipi 
satisfies ([M]) . For the grid in the left plot, we may take Vh = acpi + I3{(fj + ipk) 
for any a, /? € M, since b ■ S/vh = in the part of downwind and Vh is null 
in the rest of ^Ih- As Proposition |3] below shows, there is no nontrivial Vh & ^ 
satisfying ((64)) for the triangulation on the right-hand side of Fig. [T] 

Nevertheless, notice that as it happens with the isolated components in Fig. HI 
the cases depicted in Fig. [7] cannot occur if one designs the mesh with a strip of 
elements along F^. However, as commented at the end of Section 13.2.11 there 
are times in practice when one cannot choose the grid so that cases like those 
depicted in Fig. |4] or in Fig. [7] may occur. We now comment on the possible 
remedies to avoid nontrivial solutions of ([27l - [29)) . For example, one can con- 
nect isolated components of Vl^ by enlarging Cl^ with triangles from fi^. We 
have checked in practice that this may leave large portions of F^ as boundary 
of f2/i, and this resulted in the presence of spurious oscillations in the computed 
approximations. A better alternative is to enlarge the grid by refining those 
elements in Vij^ upwind of the isolated component. In Fig. |5]we show the result 
of dividing into four similar triangles (red or regular refinement) the two upwind 
neighbours of the isolated components in Fig. U and using longest edge bisection 
in those triangles that inherit a hanging node. We see that after the refinement 

process Clh is connected. Also, the lack of uniqueness induced by those triangles 
in Clh with a side parallel to b downwind of fi^ can be prevented by red-refining 
their upwind triangles in as it can be seen in Fig. [SI where no side parallel 
to the wind velocity b is now downwind of 17^. 

The following result states general conditions guaranteeing that ((64]) im- 
plies Vh = 0. Consider the set of points that are not downwind of 17^, that 
is, 

d{hh) ^{xetlh I {x - t6 I i > 0} n 17+ = 0}. 
Proposition 3 Let b be constant and c = 0. Then if either diClh) = Clh or for 
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Figure 8: The triangulations in Fig. |4] after regular refinement of upwind trian- 
gles closest to isolated components. The new set fi^ for b = [1, l]"^ is shadowed 
in grey and points of the new Afs sue marked with circles. 




Figure 9: The triangulations in Fig. [7] after regular refinement of triangles in 
upwind of The set fi^ for b = [1, 1]-^ is shadowed in grey and points of the 
new J\fs are marked with circles. 

any x € flh\d(flh) there is a path in Clh from x to a point y G d{Qfi) through 
elements with no edge or face parallel to b, then the only element € Vh 
satisfying ^64^ is Vh =0. 

Proof, livfi & Vh satisfies (|64|) . then —Oon d{ilh), since any Vh satisfying (|64ll 
is constant along the characteristics x + tb which, eventually intersect F" where 

Vh. = 0. If d{ilh) = dh, then Vh can only be nonzero on But elements 
on have vertices either on F^ where Vh = or on Clh- Thus, Vh must be 
zero on every element in il^. 

Suppose now that, d{(}h) ^ (ih- Then, for x e tlii\d{Uh), let 7 be the 

path in Clh connecting x to a point y e d(Clh), not intersecting any element 
with a side or face parallel to b. There is no loss of generality in assuming the 
path to be a polygonal, and except maybe from the first and last segments, the 
remaining segments joint the bary centers or arithmetic means of the vertices of 
the elements it intersects. There, wi ll be a first element r of those intersected 
by 7 which is not entirely inside Qh\d{(lh)- Since r has no edge or face parallel 

to &, then its interior f is not entirely contained in rih\d(Clh) , and thus, Vh must 
vanish on t. But then, the previous element is in a similar situation, with Vh 
vanishing on the side or face in common with r and this side or face not being 
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parallel to b, so that Vh is zero in that element too. Repeating the argument we 
conclude that Vh vanish in a;. □ 



Finally, for those situations in practice where one has to work with a mesh 
without a strip of elements along F^, we now study how to make robust the 
technique of red-refining adequate triangles in 17^. We first notice that it is 
not difficult to design examples where red-refining even all triangles in does 

not make Qh connected. However, we now argue that two red-refinements are 
enough to connect any isolated component. This is done as in the proof of Propo- 
sition [3] by considering a polygonal path in 51 joining the isolated component 
with with another component upwind of it. The segments of the path may be 
assumed to join barycenters (resp. arithmetic means of the vertices) of triangles 
with a common side (resp. tetrahedra with a common face). Consider a refine- 
ment strategy where all sides or edges are bisected, and apply it to all elements 
in 51^ intersecting the path. It is not difficult to analyze all possible cases, the 
worst cases being those where all vertices are in F^ and all sides (resp. faces) 
except those two intersecting the path are also on F^. After two refinements 
the elements in 51^ are restricted to be in the area or volume induced by the 
new vertices of the second refinement closest to the initial vertices, as we show 
in Fig. IIOI Then, in the two-dimensional case, the original path is completely 

embedded in flh, and, in the three-dimensional case, it is on the border of Cth, 

so that a small change moves it into Clh- 



Figure 10: Two and three-dimensional elements showing shadowed in grey the 
part where $7^ is restricted to be after two refinements where all sides and edges 
are bisected. Dotted lines are path entering and leaving the element through 
the only sides and faces not on the boundary of and passing through the 
barycenter. 

Let us mention that in the case of tetrahedra, the so-called regular refinement 
is not unique [3], and that the arguments above also apply if the refinement 
is done by successive bisection of the six edges of the tetrahedron using, for 
example, the technique in [1], or, in the case of triangles, three applications of 
longest-edge bisection to the three sides of a triangle. 
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4 Numerical Experiments 



In this section we solve ([THll) on different domains fl and with different forcing 
terms and vector fields b. In all cases we take c — 0. We used Matlab 7.13.0 and 
the backslash command to solve linear systems. Also, for the SUPG method, we 
take the streamline diffusion parameter as suggested in formulae (5-7) in |26| . 
which we reproduce here for the convenience of the reader. More precisely, in 
the SUPG method test functions on element r are of the form iph + Srb ■ V(/?ft, 
where S^- is given by 

^r = ^^^^^, ifPer>l, (65) 
Sr = ^^^'^('^'^) ^ if < 1^ (66) 

where Per = ''ldiam(r,b) ^ ^^^^^ </5i, . . . , V'd+i ^-re the basis functions in ele- 
ment T (taking value 1 on one of the vertices and on the rest of them) then 

diam(T, 6) — ^ ^ 



\b-Vipi\ + ■■■ + \b-Vtpd+i\' 

Here |6| stands for the euclidean norm of the vector field b. If b is not constant, 
it is evaluated at the barycenter of element r. The Matlab codes used in this 
section are available from the author on request (check also the url address 
http: //personal .us . es/bga/bgajfiles/sof tware_bga.html). Excluded are 
the codes corresponding to Example 7 which is subject of current resarch. 

Example 1 Simulation of Shishkin grids. Since the SMS method was conceived 
as a simulation of a Shishkin grid, we now check how well it does it. We take 
n = (0,1)2, ^ ^ [2,3]^, Dirichlet homo geneous boundary conditions and the 
forcing term / such that the solution is 

u{x, y)^{x- e2(--i)/-) (y2 _ e3(y-i)/^) . 

This example is taken from Example 5.2 in 25 , but in our case c — 0. On the 
one hand, we solve this problem with the SUPG method on Shishkin meshes 
with J = 2N subdivision in each coordinate direction. They are formed as 
tensor products of one-dimensional Shishkin meshes, with values ax = 2s log(7V) 
and (Ty — (3/2)elog(J) on the x and y directions, respectively. On the other 
hand, we take the coarse part of the Shishkin mesh, which is a triangulation of 

r!, = [0,1-CT,] X [0,1-aj,], (67) 

with N subdivisions in each coordinate direction, and apply the SMS method 
(both Galerkin-based and SUPG based). For iV = 5, 10, 20, . . . , 320, we compute 
the numerical approximations and measure both the computing time and the 
L°° errors in the interior points of the coarse part of the Shishkin grid (points 
which are shared by the three methods) for e = 10^^ and e = 10^^. Results 
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Figure 11: Relative efficiency of SMS and SUPG methods in Example [TJ 

are shown in Fig. [TT] Looking at the errors, and for both values of e, we notice 
that the three methods commit roughly the same errors, being in this example 
those of the SUPG-based SMS method slightly worse in general. In terms 
of computational efficiency, though, both SMS methods are roughly equally 
efficient, and markedly better than the SUPG method on Shishkin grids, since 
although the three methods commit roughly the same errors, the SMS methods 
compute the approximations between 4 and 2 times faster than the SUPG on 
the Shishkin grid (only for N = 320, the SMS methods were only 1.2 times 
faster than the SUPG method). In the next example, we further comment on 
computational costs and the structure of the linear systems to be solved to 
obtain the different approximations. 

Example 2 Comparisons on the same grid. In the previous example, it may 
be considered unfair to compare the new methods on an x iV coarse grid 
with the SUPG on a 2N x 2N grid, since this is bound to be more expensive. 
In the present example we apply the methods on the same grids (uniform with 
diagonals running Southwest- Northeast). The convection- problem is the same 
as in the previous example, and, as before, L°° errors are measured in inte- 
rior mesh points. The SUPG method in this example gave very poor results, 
with errors above 10~^. For this reason, and following suggestions in [30], we 
programmed it with some crosswind diffusion. More precisely, if in the SUPG 
method test functions on element r are of the form iph + S-rb ■ \/(ph, where dr is 
the stabilization parameter, we used (ph + Sr{b ■ Viph + ^cdxf) where the value 
of 5c was set by trial and error to obtain the smallest errors. These values 
were 5c = 0.7701, 0.8783 and 0.9365, for N = 10, 20 and 40 subdivisions in 
each coordinate direction. For similar reasons, the value of 5r in formulae (5-7) 
in |26'. was multiplied by 1.57, 1.615 and 1.64 for the above-mentioned number 
of subdivisions, respectively. 
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The results can be seen in Fig. [T^l Even though the SMS methods can be 
up to twice as costly as the SUPG on the same grid, the errors are, however, 
between 13 and 60 times smaller (SUPG-based SMS) and 26 and 130 (Galerkin- 
based SMS), so that the SMS method is computatinally much more efficient. 
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Figure 12: Relative efficiency of SMS and SUPG (with crosswind diffusion) 
methods in Example [2] 



Let us comment here on the computational cost of the methods. Recall that 
the SMS approximation Uh is obtained (together with the Lagrange multipli- 
ers Zfi and the values tj , j G as the solution of the optimality conditions 

Let ifi,. . . ,ipn be the nodal basis of Vh (each tpi takes value 1 on one single 
node of the triangulation and on the rest of them), and let A and S be the 
n X n matrices with entries 



l<i,j <n, 



(or ttij = ahi<fi,<Pj) in the case of the SUPG method) respectively, and let E 
be the n x m matrix whose columns are those of the n x n identity matrix cor- 
responding to the indexes in N5. Notice that we may assume m ~ pn^''"^^/'^ for 
some p > 0, {d being the dimension of the euclidean space where the domain 
is). Observe also that A, S, and E are typically sparse matrices. The nodal 
values of Uh, the values tj, j £ Ns and the nodal values of are then obtained 
by solving a linear system whose coefficient matrix is 



M 



S 

~E^ 
A E 



(68) 



This must compared with the SUPG method where the coefficient matrix is A, 
that is, a system of order n, whereas the SMS method is of order 2n + pn^'^^^^^'^ 
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Figure 13: The domain of Example [3] (left) and a random grid (right). 

However, as Fig. [12] shows, the errors in the SMS method are so small that 
greatly compensate for the large computational cost. Also, as Fig. [Tl] shows, 
the comparison is favourable with the Galerkin method on Shishkin meshes, 
where systems of order have to be solved. 

On the other hand notice that the change of variables Zh = —Zh, changes —A^ 
and —E'^ in M to and , respectively, so that the coefficient matrix in the 
SMS method is symmetric. This allows to use methods and software for sparse 
symmetric indefinite matrices which are generally faster than methods for gen- 
eral sparse matrices (like those in the SUPG or Galerkin methods). A study 
of the performance of the direct methods available today for sparse symmetric 
matrices can be found in [20 , where the codes MA57 12 and PARDISO, gU], 
|41] . appear as best performers. Note however that only serial versions were 
tried in [20] and that that this is an area of fast development. As for iterative 
methods, the fact that SMS approximation can be found by solving a symmet- 
ric system allows the use of three term recurrence methods like the MINRES 
method of Paige and Saunders [31] (see e. g. [S], [II], [Hj, [H], and the 
references cited therein for information on preconditioning this kind of systems). 
We remark that for the systems in the SMS method, Matlab's backslash com- 
mand seems to take no advantage of the symmetry of the systems, and, thus, 
no advantage has resulted from that symmetry in the experiments reported in 
the present paper. 

Example 3 Irregular grids on curved domains. We now consider domains 
where it may not be easy to set up a Shishkin grid. In particular, we consider 
the domain D, enclosed by the following curve: 



where a = 0.9. The curve was created by tilting 45 degrees a centered trochoid 
and altering its size in order to have the shape depicted in Fig. [131 In this 
domain we consider problem (|I}[5|) , with b — [2,3]-^ and constant forcing term 
/ = 1, with Dirichlet homogeneous boundary conditions. 

We study the behaviour of SMS methods on irregular grids, built as follows. 
For a given positive integer N, we consider Delaunay triangulations with (A^+l)^ 



lit) 



26-H7(l-sin'^(2t)) I" 1 -1 1 [ 2cos(t) 
40(2 + c^)^/2 [ 1 1 J [ 2sin(t) 



a cos(2i) 
a sin(2i) 



te[0,2n], 
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points in fi, where 4iV of them are randomly distributed on dO,. If No of these 
are on , then, we built the strip of elements along as indicated in Fig. [5j 
The remaining points are generated by first fitting a uniform grid inside and 
the outflow strip, its diameter h in the x and y direction being the value for which 
the number of points is the closest to {N — 1)^ — No, and then by displacing 
each point randomly on the x and y directions with a uniform distribution 
on [—h/3, h/3]. We show one of such grids for = 20 in Fig. [131 

For N — 40, we generated 200 random grids and on each of them we com- 
puted the SUPG and SMS approximations and their error in the convective 
derivative in tlh, that is 

Wh being each of the three approximations. The errors on the 200 grids for e = 
10~^ and £ = 10^^ are depicted in Fig. UM (marked differently for each method). 
The results have been reordered so that those of the SUPG method appear in 
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Figure 14: Errors in convective derivative on random grids for e = 10 (left) 
and 10-s (right). 

descending order. We see that the SMS methods clearly improve the errors of 
the SUPG method, specially so in the case of the SUPG-based SMS method. 
Computing the ratios of the error of the SUPG method and each of the SMS 
methods and taking the arithmetic mean, the resulting values for the Galerkin- 
based and SUPG-based SMS methods are, respectively, 9.64 and 30.19 for s = 
10~^ and 11.42 and 37.69 for e = 10"^. That is, the SMS methods commit errors 
that are, on average, between 4 and 33 times smaller than those of the SUPG 
method. For e = 10~* we repeated these computations for growing values of N, 
and the rations between the errors of the SUPG method and the SMS method 
grew with N. For example, for N = 320, these were 46.93 and 380.41 for the 
Galerkin-based and SUPG-based SMS methods, respectively. 

In Fig. [iSl where we compare the the SUPG-based SMS approximation 
that produced the largest error (left) with that of the SUPG method that pro- 
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duced the smallest error (right). We notice the typical oscillations of the SUPG 




Figure 15: The worst case of the SUPG-based SMS method (left) and the best 
case of the SUPG method (right). 

method, which are located only near the outflow boundary but are of consid- 
erable amplitude. The SMS method, on the contrary presents no oscillations. 
Similarly, in Fig. 1161 we show the Galerkin-based SMS approximations that pro- 
duced the largest and smallest errors (left and right, respectively). Both of them 
compare very favourably with the best case of the SUPG method in Fig. [T5l 
We notice however that some small amplitude oscillations can be observed in 

0.4v 0.4^ 



0.3s , . " 0.3s 




Figure 16: The worst and best cases of the Galerkin-based SMS method. 

the worst case of the Galerkin-based SMS method. We believe that these are 
due to the irregularity of the grid, since making the grids less irregular (i.e., 
smaller random perturbations) diminishes these small oscillations in the worst 
case, whereas making them more irregular increases them. It is to be remarked, 
though, that the irregularity of the grid does not affect the SUPG-based SMS 
method. Some results explaining the degradation of performance of the Galerkin 
method on irregular grids can be found in '10' and '44|. Also, further numerical 
experiments with the SMS method on irregular grids can be found in [I9j . 
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Example 4 Parabolic Layers. We solve ([T]) on 17 = (0,1)^, with e — 10~®, 
b = [1,0]"^ and / constant equal to 1 with Dirichlet homogeneous boundary 
conditions. This is a well-known test case (see e.g., [H]). The solution presents 
an exponential layer at the outflow boundary at a; = 1 and parabolic or char- 
acteristic layers along y = and y = 1 (see e. g., [3^ for a precise definition 
of these concepts). In Fig. [T7] we show the SUPG and SMS approximations 
on a 20 X 20 regular grid with Southwest-Northeast diagonals. We notice that 
whereas the SUPG approximation suffers from the typical oscillations along the 
characteristic layers, the SMS approximation is free of oscillations, both in the 
exponential layer at x = 1 and along the characteristic layers. 
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Figure 17: The SUPG (left) and the SMS approximation in Example HI 

In several techniques to reduce the oscillations of the SUPG method are 
tested. In this example, for approximations Wh computed on a regular 64 x 64 
grid, the following quantities 

osc max {u;;i(0.5, y) — Wh(0.5, 0.5)}, (69) 

ye{l/64,... ,63/64} 

smear := max {u;/i(G.5, 0.5) — Wh(0.5, y)}, (70) 

iyG{l/64,... ,63/64} 

are computed in [26] as a measure of the oscillations and the smearing along the 
characteristic layers, desirable values being, respectively, between and 10~^ 
and between and 10~^. The value of osc and smear for all the methods tested 
in [35] except one was always larger than 10^*^. The values in the case of SMS 
methods were below 10"". The value of osc in the SUPG method in our tests 
coincided with that in [26 , 0.134 (no value of smear was given in ^26;). 

Similar striking contrast between methods tested in [26] and the SMS meth- 
ods can be found in the experiments on randomly-generated grids for this ex- 
ample in 19 . 

Example 5 Interior layers. This example is also taken from [26] . We solve ^ 
on f7 = (0,1)2, .^^ii-j^ ^ ^ ;^Q-8^ 5 ^ [cos(-7r/3),sin(-7r/3)]^, / = 0, and u = y 
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on dQ where 



0, ifx = lory<0.7 

1, otherwise. 



The solution possesses an interior layer starting at x = and y = 0.7, and an 
exponential layer on a; = 1 and on the right part of y = 0. For SMS methods, 
the only way we have conceived so far to deal with interior layers is to treat 
them as parabolic layers. For this to be possible, the grid has to include the 
characteristic curve that starts at the discontinuity of the boundary data (or a 
polygonal approximation to it in case of curved characteristics) . We will refer to 
this characteristic curve as the layer characteristic. For grids as those depicted 
on the left of Fig. [HI where the layer characteristic is not part of the grid, one 




Figure 18: Left, a 8 x 8 uniform grid in Example [5l with the characteristic curve 
of the internal layer (dashed line). Right, the grid on the left but moving to the 
internal layer characteristic its closest points on the grid. 

possibility to make it part of it is to move to the layer characteristic its closest 
points in the grid, as we show on the right plot in Fig. [TH] (see also the next 
example for an alternative). Once this is done, on has to set the value of the 
solution along the layer characteristic as part of Dirichlet boundary conditions. 
To do this, we integrate the reduced problem (i.e., when e = 0) along the 

layer characteristic. In the present case, this procedure gives u = 0.5 on the 
layer characteristic. 

With these provisions, we compute the approximations of the SUPG method 
SMS methods on a 16 x 16 grid. The results are shown in Fig. [19] (the SMS 
approximations were identical and only the Galerkin-based one is shown). As 
before, the SMS method produces an approximation with no oscillations, in 
sharp contrast with the SUPG method. Comparison with the methods tested 
in [26] can be found in [19], with results very similar to those of the previous 
example. 

Example 6 Hemker problem. Here, il = (-3,9) x (-3,3)\{a;2 + < 1}, 
b — [1, O]-'" and / = 0. The boundary conditions are 



u{x,y) 




1, 

eS/u • n = 0, 



if a; = 

if + = 1, 
elsewhere. 



~3, 

,2 - 



(71) 
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The SUPG approximation 



Galerkin-based SMS approximation 




Figure 19: The SUPG and SMS (Galerkin-based) approximations on a 20 x 20 
regular grid in Example O 

This problem, which was originally proposed in [22] , models a hot column (the 
unit circle) with normalized temperature equal to 1, and the heat being trans- 
ported in the direction of the wind velocity b. Thus, a boundary layer appears 
in the upwind part of the unit circumference from the lowest to highest point, 
and two internal layers start from these two points and spread in the direction 
of b. Notice also that part of the boundary is curved, a feature which is often 
encountered in applications. 

We are going to present results corresponding to e = 10^^ on the grid 
shown in Fig. [20l which has 932 triangles and 531 nodes. The interior layer 




-2 2 4 6 8 



Figure 20: The grid in Example [6) 

characteristics are not part of this grid, as it can be seen in Fig. [3TJ where 
we show one of the layer characteristic in a dashed line. We have seen in the 
previous example that interior layers are treated in the SMS method in the same 
way as characteristic boundary layers, so that they must be part of the grid. 
Rather than moving grid points to the layer as we did in the previous example, 
in the present one we will enlarge the grid with more triangles and vertices so 
that the layer characteristic is part of it (an example is shown in Fig. This 
may be useful in those cases in practice where, as commented at the end of 
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Section [3.2.1l one does not have the freedom to choose or move the grid points. 

Thus, we now describe a general technique to enlarge grids in order to include 
an inner layer characteristic. For simplicity we describe it for two-dimensional 
problems. Let 7 be the inner layer characteristic. We assume that the mesh is 
sufficiently fine so that 7 can be well approximated by a straight segment in the 
interior of every triangle r it intersects. Let r be such a triangle. If 7 passes 
through a vertex v, then the r is bisected by 7 into two triangles. These two 
triangles are included in the enlarged triangulation. Otherwise, 7 intersects two 
sides, and the element r is divided by 7 into a triangle and a quadrilateral (see 
an example in Fig. I22p which in turn is divided into four triangles. We remark 




Figure 21: Left. Detail on the grid in Fig. [20l showing one of the inner layer 
characteristic (dashed line). Right: the same grid enlarged with more triangles 
and nodes to include the layer characteristic. New sides are plotted with thinner 
lines. 




Figure 22: Detail of the enlarged grid in Fig. [?T] showing a triangle of the original 
grid being divided by the layer characteristic in a triangle and a quadrilateral, 
which we divide into triangles by joining the vertices with their arithmetic mean. 

that it is easy to conceive better strategies to enlarge the triangulation in order 
to include the layer characteristic (techniques that avoid long-shaped triangles, 
for example), but we have chosen this one due to its simplicity. Nevertheless, as 
the experiments below show, its simplicity does not prevent it from obtaining 
excellent results. 

Fig. [23l shows the SUPG and the SMS approximations (Galerkin-based). Os- 
cillations can be clearly seen in the SUPG approximation, but they are absent in 
the SMS approximation. This can be better observed in the bottom plots, where 
a different point of view is taken. Following [2] , we measure the undershoots 
of an approximation Wh as min{u;/i} and the overshoots as max{w/i — 1}. The 
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Figure 23: The SUPG approximation (left) and the Galerkin-based SMS ap- 
proximation (right) on the grid depicted in Fig l20l A different point of view is 
taken on the bottom plots. 

over and undershoots for the SUPG method were 0.04 and —0.52, whereas for 
both SMS methods the overshoots where smaller than 10~^^, and the under- 
shoots were 0. Similar results were obtained by doubling and multiplying by 
four the number of subdivisions in each coordinate direction. The results of the 
SMS method compare very favourably not only with the SUPG method, but 
with most of the methods tested in , which had values very similar to those 
of the SUPG method. They also compare very favourably with results in [18], 
where LPS methods present oscillations of about 5% of the jump, whereas the 
SMS methods of less than 10~^^%. 

In [19j . the previous experiment is repeated with the vector field b changed 
to 6 = [cos(6'), sin(0)]^, for 100 equidistant values of 9 between (0,7r/4], with 
results similar to those shown here for all cases except four, where, in order to 
get over and undershoots of order 10~^^ it was also necessary for nodes 0(/i^ji,j) 
away from interior layers to be allowed to move (see [TH] for details). 

Example 7 Double-glazing test problem [13j. This is an example where the 
vector field b has vortices. The domain is = (—1, 1)^, / = 0, and the wind 
velocity is [y{l — a;^), —x{l — y^)] , so that the characteristic curves are the closed 
curves given by 

{{x,y) I (1 - x^)(l - y^) = constant}, 

(See Fig. [24)1 . The Dirichlet boundary conditions are u = 1 on a; = 1 and 
u — otherwise, so that there are discontinuities in the data at corners in 
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X — 1, y = ±1. These discontinuities give raise to parabolic boundary layers 
which, according to |13j . have a structure difficult to compute by asymptotic 
techniques. 




Figure 24: Left: The streamlines in Example [71 Center: the solution in Exam- 
ple [7] for e = 0.005. Right: A regular mesh, with the set Xl^ shadowed in grey 
and points of Afs are marked with circles. 

Notice that the hypotheses stated at the beginning of Section 12.21 {b{x) ^ 
and all characteristics entering and leaving the domain in finite time) do not 
apply to vector fields with vortices or closed integral curves. Nevertheless, as we 
now show, the results of the SMS method in the present example are as good 
as in the previous ones, even though part of the analysis in Section [3. 21 does not 
apply to the present case. 

Observe that since the four sides of the boundary are themselves character- 
istic curves, we have — Sfi, so that building as described in Section [2?2l 
results in $7^ being composed by all elements touching the boundary. In this 
case the SMS method produces an approximation equal to on all nodes except 
on those on a; = 1, where it takes value one, and this is only correct for s = 0. 
For £ > 0, better results are obtained with the SMS method if, as we show 
in Fig. [24l the set fi^ is shrunk so as to correspond to that of a slightly smaller 
rectangular domain, that is, elements are included in il'^ if they intersect the 
outflow boundary of (—1 + 5,1 — 5)^ for a small 6 > (other possibilities to 
obtain better results with vector fields with vortices will be reported in fu- 
ture works). With this selection of fi^, the results of the Galerkin-based SMS 
method for e = 10~'^ can be seen in Fig.[25l Also shown in Fig. [25]is the SUPG 
approximation, where we can observe oscillations along the characteristic layers, 
especially at a: = 1 and y = —1. In sharp contrast, the SMS method produces a 
nonnegative approximation. For this value of e, in order to obtain nonnegative 
approximations on a regular N x N mesh with the SUPG method one must take 
N > 144, and N > 240 with the Galerkin method; for e — 10^^, this is achieved 
for N > 682 and N > 800 with the SUPG and Galerkin methods respectively, 
and for e = 10~^, neither method was we able to obtain nonnegative approxi- 
mations with all the meshes we tried up to iV = 1800. In sharp contrast, with 
the SMS method, both Galerkin-based and SUPG-based, nonnegative approxi- 
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Galerkin-based SMS SUPG 




Figure 25: The Galerkin based SMS method (left) and the SUPG method on a 
20 X 20 regular grid in Example [7] for e = 10^"'. 

mations were obtained for > 8 and e as small as 10"^''. 

5 Concluding remarks 

A novel stabilization technique for convection-diffusion problems has been intro- 
duced, tested and partially analyzed. It can be applied to most existing methods 
based on conforming piecewise linear finite elements. It consists of, first, adding 
extra values to the residual of the method on nodes next to the outflow and 
characteristic boundaries. The resulting equations are taken as a restriction on 
a least-squares problem on the elementwise residuals of the convection-diffusion 
operator, where elements next to the outflow and characteristic boundaries (and 
internal layers) are left out. 

The method has been tried in a series of standard and nonstandard tests, and 
results seem to suggests it performs manifestly better than a good deal of the 
methods of choice today. The tests include exponential and characteristic layers, 
and even irregular grids on domains with nontrivial geometries. This is so in 
spite of the method being initially conceived as a simulation on the coarse part 
of a Shishkin mesh. The tests also include interior layers and convection with 
vortices, with similarly outstanding results, although numerical experiments in 
the present paper and in jl9j suggest that further research may be needed on 
these topics. 

Besides the practical performance shown in tests, it is to be remarked the 
lack of parameters in the method (in markedly contrast with most of existing 
stabilized methods today). Subject of future research will be extending the 
method to finite elements of higher degree, as well as the possible changes when 
the mesh Peclet number tends to one. 

Acknov^rledgements. The author wish to thank Prof. Martin Stynes for advice 
and helpful discussions in the research summarized in this paper. 
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Shishkin mesh simulation: Complementary 
experiments 

Bosco Garcia- ArchillaH 
January 18, 2013 

Abstract 

We present here some of the experiments that, for the sake of brevity 
were not included [3]. The present paper is neither self-contained nor 
intended to he published in any major journal, but is intended to be ac- 
cessible to anyone wishing to learn more on the performance of the SMS 
method. It should be read as a complement to [3]. 

1 Numerical experiments 

In its present version, this document is not selfcontained, and its contents must 
be seen as add ons to [3]. Recall that in that paper we are concerned with 

- eAu + 6 • Vu + cu = /, in fi, (1) 

dii 

u = gi, indflD, t;- = 32, in 917 at, (2) 
on 

where f7 is a bounded domain in M'*, d = 1,2,3, its boundary dft being the 
disjoint union of F/) and Tn, b and c are given functions and e > is a constant 
diffusion coefficient. 

Example 1 Simulation of Shishkin grids. Let us also mention that we rerun 
the SUPG-bascd SMS manually tunning the streamline diffusion parameter and 
found that, in this example, if set 1.5 times larger than the value shown in |31 
Fig. 11], the errors improve by a factor of two approximately, resulting in the 
most efficient method and at least twice as efficient as the SUPG on Shishkin 
grids 

Example 2 Comparisons on the same grid. The difference between the SUPG 
and the SMS methods is even larger when errors are meassued in the norm, 
which are shown in Fig. [T] 

Example 3 Irregular grids on curved domains. 

■^Departamento de Matematica Aplicada II, Universidad de Sevilla, Sevilla, Spain. Re- 
search supported by Spanish MEC under grant MTM2009-07849 {bosco@esi.us.es) 
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Figure 1: Relative efficiency of SMS and SUPG (with crosswind diffusion) metli- 
ods in tfie /f^norm. 



We study the behaviour of SMS methods on irregular grids. For a given 
positive integer N, we consider grids with {N + 1)^ points in fl, where AN of 
them are randomly distributed on dfl. Let us denote by No the number of these 
on r^. Then, in from the remaining (TV — 1)^ points, No of them are displaced 
by the normal vector from those on (this is done to create a strip of elements 
on the outflow boundary) and the rest are generated randomly according to the 
following two schemes. 

Highly irregular grids The points are generated following uniform distribution 
on the square enclosing ft, and those enclosed by dil and the outflow strip 
are included in the grid until the desired number {N — 1)^ — No of them 
is reached. 

Mildly irregular grids A uniform grid is build flrst inside fl and the outflow 
strip, its diameter h in the x and y direction being the value for which 
the number of points is the closest to {N — 1)^ — Nq. Then each point is 
displaced randomly on the x and y directions with a uniform distribution 
on [—h/3, h/3\. This is the type of grid used in the experiments in [3]. 

An example with the two types of grids can be seen in Fig. |4l 

For different values of N, we generated 200 random grids of both types, and 
on each of them we computed the SUPG and SMS approximations and their 
error in the convective derivative in 17^, that is 



\b-\/wh - 11 



for Wh = Uh the SUPG approximation and wt = u/i, the SMS approximation 
(both Galerkin and SUPG-based). The errors on the 200 grids for iV = 80 
and e = 10~^ arc depicted in Fig. [2] They are marked with a dot, an asterisks 
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and a cross for the SUPG, Galerking-based SMS and SUPG-based SMS meth- 
ods, respectively. The results have been reordered so that those of the SUPG 
method appear in descending order. As in the results shown in [3], the SMS 
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Figure 2: Errors in convective derivative on random grids. 

methods clearly improve the error with respect the SUPG method, particularly 
so in the case of the SUPG-based SMS method. Computing the ratios of the 
error of the SUPG method and each of the SMS methods and taking the arith- 
metic mean, the resulting values are for the Galerkin-based and SUPG-based 
SMS methods, respectively, 5.17 and 65.57 on highly irregular grids, and 16.30 
and 75.16 on mildly irregular grids. That is, the SMS methods commit errors 
that are, on average, between 5 and 75 times smaller than those of the SUPG 
method. 

In Fig. [5] we can also see that while the SUPG and SUPG-based SMS meth- 
ods are insensitive to the irregularity of the grids, this is not the case of the 
Galerkin-based SMS method. This is also reflected on Fig. [3l where we show 
the worst cases of Galerkin-based SMS method on both types of random grids 
for = 40: whereas in the highly irregular grids the approximation exhibits 
oscillations, this not so noticeable on the mildly irregular grid (the correspond- 
ing grids are shown in Fig.^l). As commented in 131 , some results explaining the 
degradation of performance of the Galerkin method on irregular grids can be 
found in [Tj and [5]. However, we must remark that the approximation on the 
left plot in Fig. |3]is the worst case out of 200, and that in most of the cases on 
highly-irregular grids the Galerkin-based SMS method did not present spurious 
oscillations. Nevertheless, due to the better performance of the Galerkin-based 
SMS method on mildly-irregular grids, in the rest of the paper, we will only 
deal with cither regular or mildly-irregular grids. 

We show the arithmetic mean of the errors committed by the different meth- 
ods on mildly irregular grids for each value of A^ in Fig. [5] We also show the 
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Figure 3: The worst cases of the Galerkin-based SMS method for A'' = 40 on 
highly irregular and mildly irregular grids (left and right, respectively). 

slopes of least-squares fits to the different sets of results. We observe that the 
errors in the SUPG-based SMS method decay as the meshes become finer (at a 
rate of N~^-^) marginally so in the case of the Galerkin-based SMS method and 
they grow in the SUPG method. That is, the SUPG method has not reached 
the convergence regime yet for these meshes in this problem. We can have an 
idea of why the errors grow by looking at Fig. [51 where we show the SUPG 
approximation on the grids shown in Fig. 2] (compare with the Galerkin-based 
SMS in Fig. [S]). We notice that large amplitude of the oscillations near the 
outflow boundary. The large amplitude of the oscillations decreases slower than 
the grid size, so that the errors in the convective derivative grow. 

Example 4 Parabolic Layers. Recall that in [3] we solve ((TJ on = (0, 1)^, 
with e = 10^*, b ~ [1, 0]-^ and / constant equal to f with Dirichlet homogeneous 
boundary conditions. The solution has one exponential layer at a; = 1 and two 
parabolic layers along ?; = and y = 1- 

In [3], a random grid is used and the sets ^2 = (0, 0.9) x (0, 0.1] and fl^ = 
(0, 0.9) X [0.9, 1], and ^4 = [0.9, 1) x (0.1, 0.9) are considered, and the following 
quantities computed for numerical approximations Wh' 

0SCpara(2) := max{ max -dyWh[xs,ys), max dyWh{xs,ys)] , (3) 

osccxp := max dxWh{xs,ys) (4) 

(£Cs, 1/3)604 

{xs, ys) being the barycenters of the triangles. The optimal values for 0SCpara(2) 
and osCoxp are, respectively, and 1, the larger these values are, the stronger are 
the oscillations in the characteristic and exponential layers, respectively. None 
of the methods tested in [3] presents a value of 0SCpara(2) below 10^^, nor a 
value of osCexp that is less than 10~^ away from the target value 1. We try 
200 random irregular grids with 1681 interior grid points (the grid in |5) had 
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Figure 4: The grids of the plots in Fig. |3l The set is shadowed in grey and 
the points in A/5 are marked with a dot. 

1721), which were created as in Example |3l On these grids, the highest values 
of 0SCpara(2) ^nd |osCcxp — 1| in the SMS methods were, respectively, 2.9 x 10"^'^ 
and 2.8 x 10^^'^. Similar striking differences were observed in the rest of the 
quantities measured in in this example. 

Also, with highly irregular grids, as long as they had a strip of elements 
along r'^, the results were similar to those of mildly- irregular grids. 

Example 5 Interior layers. Nothing to be added here at this stage. This 
example is revisited after Exameple IHl below. 

Example 6 Hemker problem. We test the technique to include the layer char- 
acteristic in the grid that was proposed in '3, Example 6]. To do this, we change 
b and the boundary conditions to 



for e (0,7r/4]. (The boundary conditions are changed so that the inflow 
boundary is part of the Diriclilet boundary To)- The idea is to try different 
combinations of characteristic layers and grid alignments. In the following com- 
putations, the enlarged grids and the corresponding sets £7^ can be as irregular 
as those depicted in Fig. [T] To prevent elements with one side parallel to b 
being downwind of fi^ (see [21 Section 3.2.2] for the possible adverse effects of 
this) when 9 = ir/A we red-refined all triangles in along the interior layers. 
Otherwise, no extra provisions where taken in the computations that follow, 

except removing isolated components in fl^ (when present) as explained in [3l 
Section 3.2.2]. 



b = 



cos{e) 
sm{e) 




if X — —3 or y 
if x^ + y'^ = 1, 



elsewhere 
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Figure 5: Arithemetic means of the errors on mildly irregular grids. 




Figure 6: The SUPG approximation on the grids shown in Fig. |4l 

We considered 100 equidistant values of 6 in (0,7r/4]. For each of them we 
computed the overshoots and undershoots of the SUPG and SMS approxima- 
tions on the grid depicted in [3l Fig. 20]. For the SUPG and the Galerkin-based 
SMS approximations, Fig. |8] shows the difference between the overshoots and 
undershoots, as a measure of the oscillations in each method (the results of the 
SUPG-based SMS method were similar to those of the Galerkin-based method 
and are not shown). We see that out of the 100 values of 9 tried, only in four of 
them are the oscillation on the SMS method larger than 10~^^, the arithmetic 
mean of these four values being 0.013, yet well below the average value in the 
SUPG method, 0.65. 

We investigate the reason of those four cases presenting values of over and 
undershoots so different to the rest of the cases. The reason is that some mesh 
points of the enlarged grid are too close to the interior layer, and this results 



6 



Figure 7: Detail of the enlarged grid, the set fi^ (shadowed in grey) and the 
points in Ms (marked with asterisks) in Example [5] for = 7r/6 in ([S]). 
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Figure 8: Difference between the over and undershoots for the approximations 
on the grid depicted in [3, Fig. 20] in Example|6l with h and boundary conditions 
as specified in ([5]). 

in the set fi^ being too thin in the neighbourhood of those points. This can 
be seen in Fig. |9l where on the left plot we show a detail of the enlarged grid 
(with shadowed in grey) for 9 — 237r/400, which is the first value in Fig. [5] 
where the overshoots or undershoots are above 10^^^ in the SMS method. The 
center plot is a magnification of the left one around the point too close to the 
layer characteristic, where we can see that very stretched triangles have been 
created on the enlarged grid. If we move the point to the layer characteristic, 
the resulting set is shown in the right plot of Fig. [SI and we can see it is 
much wider than that in the left plot. 

We repeated the four cases in the SMS method were oscillations were much 
larger than the rest of the cases, but this time allowing a vertex v oi a. triangle r 
to move to the layer characteristic if its distance to it was less than (/imin(T))^/10, 
were hniin{T) is the minimum distance between the vertices of r. Once the SMS 
approximations were computed on the moved grids, values on the unmoved grid 
were computed by interpolation and over and undershoots were again computed. 
The overshoots in all cased turned out to be below 10"^*, and the undershoots 
were zero in all cases. Notice that when not altering the grid (not even by less 
that /i^jjj/lO) is an important issue, then, another option to avoid mesh points 
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Figure 9: Left: detail of the enlarged grid and for 9 — 237r/400 in Example[S] 
with ([5]). Center: magnfication of the left plot around a point too close to the 
characteristic layer. Right: the set after moving that point to the layer 
characteristic. 

too close to the layer characteristic is to slightly perturb the wind velocity b in 
the neighbourhood of those points so that the layer characteristic passes through 
them. 

Example 5 (revisited). As in example HI we compute with our method some 
of the quantities used in [3] to measure the quality of the approximations. For 
an approximation Wh € Vh, these are 



OSCint : = 

smearjut := 




where fii — {{x,y) G ft \ x < 0.5, y > 0.1}, and xi and X2 are the 
first points on the line y = 0.25 satisfying, respectively, Wh{xi,0.25) > 0.1 
and 'Wh{x2,0.25) > 0.9. It is argued in [4] that ([7]) is a measure of the thickness 
of the interior layer. Observe also that the target value for osCint is 0. 

We computed the value of these quantities in the case of the SMS approxi- 
mations on a regular 64 x 64 grid with Southwest-Northeast diagonals, using the 
technique of the previous example to include the layer characteristic in the grid. 
On a similar mesh, out of 18 methods methods tested in [J] the value of osCint 
was larger than 0.1 in 14 of them, and larger than 0.001 in three more of them. 
In contrast, the value of oscjnt in the SMS methods was 1.8 x 10"^'* (Galerkin- 
based) and 2.2 x 10^^'^ (SUPG-based). Also, the vale of smearjnt was larger 
than 6.2 x 10~^ in all the methods tested in whereas for both SMS methods 
was 1.3 X 10~^. For the SUPG method, the values of osCint and sniearint were, 
respectively, 0.59 and 0.062, the same values as in [4j. 

On the irregular grid in 4 , they obtain results were very similar to those of 
the regular grid. We, on our part, tried the 200 random grids of Example 4 in 
the present one. The mean value smearint in the SMS methods was 2.7 x lO^'^. 
With respect the value of osCint in the SMS methods, out of the 200 runs, in only 
20 of them was this value larger than 10^*, being the mean value in these 20 
cases of 0.1822, which still compares very favourable with most of the methods 
tested in |4] . We rerun these 20 cases allowing a grid points to move to the layer 



8 



characteristic if the distance to it was less 'ih'^-^, and in 19 of them the value 
of oscint move to under 10~^°. In the remaining case this was also achieve if we 
red-refining the triangles in $7^ next to the inner layer characteristic. 

Example 7 Double-glazing test problem Whatever we may have to add in 
this case will be reported elsewhere. 
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